rm(list = ls())
library(countrycode)
library(pscl)
library(foreign)
library(readstata13)
library(wnominate)
library(reshape)
library(polycor)
library(ggplot2)
library(ggthemes)
library(corrgram)

setwd("~/Desktop/QJPS Replication")

for(iii in 1:13){
  
  data <- read.csv("data/roll_call_matrix.csv")
  
  ##Miller et all democracy data
  dem <- read.dta13("data/ccpcnc_v2_small.dta")
  
  #combine democracy data with roll call matrix
  data$match <- paste(data$cowcode, data$year)
  dem$match <- paste(dem$cowcode, dem$year)
  m <- match(data$match, dem$match)
  table(is.na(m))
  
  dem$country <- as.character(dem$country)
  data$country_name <- dem$country[m]
  data$country_name <- as.character(data$country_name)
  
  #get rid of Mexico
  data <- data[data$country_name != "Mexico",]
  
  #subset data to only those that had a change occur
  data2 <- data[data$evnt == 1,]
  data2 <- data2[,11:(ncol(data2)-2)]
  
  #How many missing values are there for each entry
  num.missing <- apply(data2, 1, function(x)sum(ifelse(is.na(x), 1, 0)))
  prop.missing <- num.missing / ncol(data2)
  # hist(prop.missing)
  
  #Get rid of those entries with more than 90 percent missing
  data3 <- data2[prop.missing < .9,]
  
  group <- read.csv("data/groups_by_vote_type.csv")
  gp <- group$group
  
  #omit a group category and make the roll call matrix
  data3 <- data3[,gp != iii]
  
  #make the matrix a roll call object
  mat <- rollcall(data3, yea = 1, nay = 0, legis.names = rownames(data3))
  
  ##########################
  #scale that beautiful roll call matrix
  set.seed(4227492)
  idealpoints <- wnominate(mat, dims = 1, polarity = 1, lop = 0.05)
  
  points <- as.data.frame(idealpoints$legislators$coord1D)
  
  names <- as.character(data$country_name[data$evnt == 1][prop.missing < .9])
  ccode <- as.character(data$cowcode[data$evnt == 1][prop.missing < .9])
  year <- as.character(data$year[data$evnt == 1][prop.missing < .9])
  
  points$name <- names
  points$ccode <- ccode
  points$year <- year
  names(points) <- c("coord1D", "Name", "cowcode", "year")
  
  ########SAVE THE DATA
  write.csv(points, paste("omitted_", iii, ".csv", sep = ""))
  
  print(iii)
}




data <- read.csv("all_countries_ideal.csv")
data$code.year <- paste(data$cown, data$year, sep = "-")

data1 <- read.csv("omitted_1.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit1 <- data1$coord1D[m]

data1 <- read.csv("omitted_2.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit2 <- data1$coord1D[m]

data1 <- read.csv("omitted_3.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit3 <- data1$coord1D[m]

data1 <- read.csv("omitted_4.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit4 <- data1$coord1D[m]

data1 <- read.csv("omitted_5.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit5 <- data1$coord1D[m]

data1 <- read.csv("omitted_6.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit6 <- data1$coord1D[m]

data1 <- read.csv("omitted_7.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit7 <- data1$coord1D[m]

data1 <- read.csv("omitted_8.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit8 <- data1$coord1D[m]

data1 <- read.csv("omitted_9.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit9 <- data1$coord1D[m]

data1 <- read.csv("omitted_10.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit10 <- data1$coord1D[m]

data1 <- read.csv("omitted_11.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit11 <- data1$coord1D[m]

data1 <- read.csv("omitted_12.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit12 <- data1$coord1D[m]

data1 <- read.csv("omitted_13.csv")
data1$code.year <- paste(data1$cowcode, data1$year, sep = "-")
m <- match(data$code.year, data1$code.year)
table(is.na(m))
data$omit13 <- data1$coord1D[m]

###############
## Correlation Matrices

rc.matrix <- data[,c(2, 10:ncol(data) )]

corr.matrix <- rc.matrix
names(corr.matrix) <- c("all features", "general", "amending", "executive", "legislature", "judiciary", "federalism", "elections", "regulatory", "international", "duties", "criminal", "rights", "special")

cor.matrix <- round(cor(corr.matrix, use = "complete.obs", method = "spearman"), digits = 3)
cor.matrix[2,1] <- NA
cor.matrix[3,1:2] <- NA
cor.matrix[4,1:3] <- NA
cor.matrix[5,1:4] <- NA
cor.matrix[6,1:5] <- NA
cor.matrix[7,1:6] <- NA
cor.matrix[8,1:7] <- NA
cor.matrix[9,1:8] <- NA
cor.matrix[10,1:9] <- NA
cor.matrix[11,1:10] <- NA
cor.matrix[12,1:11] <- NA
cor.matrix[13,1:12] <- NA
cor.matrix[14,1:13] <- NA

cor.matrix <- abs(cor.matrix)
cor.matrix

cor.dat <- melt(cor.matrix)
cor.dat <- cor.dat[-which(is.na(cor.dat[, 3])),]
cor.dat <- data.frame(cor.dat)
cor.dat

levels(cor.dat$X1) <- list("all features" = "af", "general" = "gen", "amending" = "amnd", "executive" = "exec", "legislature" = "leg", "judiciary" = "jud", "federalism" = "fed", "elections" = "elec", "regulatory" = "reg", "international" = "int", "duties" = "dut", "criminal" = "crim", "rights" = "rgt", "special" = "spec")
levels(cor.dat$X2) <- rev(list("all features" = "af", "general" = "gen", "amending" = "amnd", "executive" = "exec", "legislature" = "leg", "judiciary" = "jud", "federalism" = "fed", "elections" = "elec", "regulatory" = "reg", "international" = "int", "duties" = "dut", "criminal" = "crim", "rights" = "rgt", "special" = "spec"))

theme_set(theme_bw())

ggplot(cor.dat, aes(X2, X1, fill = value)) + 
  geom_tile() + ggtitle("Correlations of Models Run with Omitted Items") +
  geom_text(aes(X2, X1, label = value), color = "black", size = 2) +
  scale_fill_gradient(name=expression("Spearman" * ~ rho), low = "#ffffff", high = "#045975", breaks=seq(0, 1, by = 0.2), limits = c(0, 1)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  labs(x = "", y = "") + 
  guides(fill = guide_colorbar(barwidth = 10, barheight = 1, title.position = "top", title.hjust = 0.5)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), 
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(1, 0),
        legend.position = c(0.9, 0.7),
        legend.direction = "horizontal") +
  guides(fill = guide_colorbar(barwidth = 7, barheight = 1, title.position = "top", 
                               title.hjust = 0.5))

mean(cor.dat$value[cor.dat$value != 1])



